Mise en place

Chargement des packages:

library(lme4) # pour la modélisation multiniveau
library(multilevelTools) # pour l'analyse des modèles
library(lmerTest) # pour l'évaluation des modèles
library(sf) # pour les objets spatiaux (cartes)
library(tidyverse) # pour la manipulation des données
library(readxl) # pour ouvrir les fichiers xlsx
library(ggrepel) # pour des étiquettes de graphiques propres

Chargement des données:

Jeu de données de niveau 1

Il s’agit des vignobles de Bourgogne. On a dans un tableau les caractéristiques de 2391 vignobles en termes de:

prix moyen du vin (millésimes 1998, 2002 et 2003)
qualité pédologique et météo (radiations solaires, pluie, etc.)
appellation d’origine contrôlée en 5 niveaux (Côteaux bourgignons < Bourgogne < Village < Premier cru < Grand cru.)
surface

Sources: Ay J.S., Hilal M., 2020, « Les déterminants naturels et politiques des AOC viticoles de Côte-d’Or », Cybergeo: European Journal of Geography, document 973, DOI : https://doi.org/10.4000/cybergeo.36443 et https://www.hachette-vins.com/guide-vins/actualite-vin/487/les-meilleurs-millesimes-des-vins-de-france/

vignobles <- st_read("Data/Geo/Vignobles_prix.shp",
                     stringsAsFactors = T) %>% 
  st_transform(crs = "EPSG:2154")

## Reading layer `Vignobles_prix' from data source 
##   `/Users/ccottineau/Documents/GitHub/QuantiLilleGeo/Cours_VinBourgogne/Data/Geo/Vignobles_prix.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2391 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 826631.6 ymin: 6645988 xmax: 851465.7 ymax: 6690810
## Projected CRS: RGF93 v1 / Lambert-93

summary(vignobles)

##                              Concat                     LIBCOM    
##  1AUXEY-DURESSESCORLOIN         :   1   MEURSAULT          : 158  
##  1AUXEY-DURESSESDERRIERE LE BOIS:   1   GEVREY-CHAMBERTIN  : 128  
##  1AUXEY-DURESSESEN CHAIGNUT     :   1   BEAUNE             : 126  
##  1AUXEY-DURESSESEN FRIVOLE      :   1   NUITS-SAINT-GEORGES: 121  
##  1AUXEY-DURESSESEN MORVIN       :   1   SAINT-AUBIN        : 121  
##  1AUXEY-DURESSESEN POCHENOT     :   1   (Other)            :1727  
##  (Other)                        :2385   NA's               :  10  
##             NOM               NIVEAU        SURFACE          SCORE_b     
##  LE VILLAGE   :  26   Bourgogne  : 624   Min.   :0.0000   Min.   :18.40  
##  LES CRAIS    :  10   Coteaux b. : 319   1st Qu.:0.1300   1st Qu.:67.90  
##  LES CRAS     :   8   Grand cru  :  37   Median :0.3100   Median :73.20  
##  SOUS LA VELLE:   8   Premier cru: 410   Mean   :0.4773   Mean   :71.89  
##  LES PERRIERES:   7   Village    :1001   3rd Qu.:0.6100   3rd Qu.:76.60  
##  (Other)      :2331                      Max.   :7.2000   Max.   :94.22  
##  NA's         :   1                                                      
##             Source        PrixMoy                 geometry   
##  guide_hachette: 416   Min.   :    5.0   MULTIPOLYGON :2391  
##  interpolation :1942   1st Qu.:   23.0   epsg:2154    :   0  
##  wine_search   :  33   Median :   36.0   +proj=lcc ...:   0  
##                        Mean   :  154.4                       
##                        3rd Qu.:   87.0                       
##                        Max.   :25960.0                       
##                        NA's   :4

vignobles %>% 
  dplyr::filter(LIBCOM == "FLAGEY-ECHEZEAUX" & NIVEAU == "Premier cru")

Jeu de données de niveau 2

Il s’agit des communes de Côte d’Or. On a dans un tableau les caractéristiques de 31 communes en termes de:

superficie
population
côte d’appellation
hierarchie administrative

Source: https://geo.data.gouv.fr/en/datasets/cac9f2c0de2d3a0209af2080854b6f6a7ee3d9f4

communes <- st_read("Data/Geo/communes.shp",
                     stringsAsFactors = T) %>% 
  st_transform(crs = "EPSG:2154")

## Reading layer `Communes' from data source 
##   `/Users/ccottineau/Documents/GitHub/QuantiLilleGeo/Cours_VinBourgogne/Data/Geo/Communes.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 31 features and 20 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 826428 ymin: 6645899 xmax: 855115 ymax: 6691133
## Projected CRS: RGF93 v1 / Lambert-93

summary(communes)

##       gid          id_geofla       code_comm    insee_com 
##  Min.   :  563   Min.   :  563   010    : 1   21010  : 1  
##  1st Qu.: 6744   1st Qu.: 6744   037    : 1   21037  : 1  
##  Median :14409   Median :14409   054    : 1   21054  : 1  
##  Mean   :11410   Mean   :11410   110    : 1   21110  : 1  
##  3rd Qu.:14602   3rd Qu.:14602   133    : 1   21133  : 1  
##  Max.   :14749   Max.   :14749   150    : 1   21150  : 1  
##                                  (Other):25   (Other):25  
##                  nom_comm               statut     x_chf_lieu     y_chf_lieu   
##  ALOXE-CORTON        : 1   Chef-lieu canton: 3   Min.   :8291   Min.   :66472  
##  AUXEY-DURESSES      : 1   Commune simple  :27   1st Qu.:8349   1st Qu.:66570  
##  BEAUNE              : 1   Sous-préfecture : 1   Median :8431   Median :66660  
##  BROCHON             : 1                         Mean   :8419   Mean   :66676  
##  CHAMBOLLE-MUSIGNY   : 1                         3rd Qu.:8486   3rd Qu.:66774  
##  CHASSAGNE-MONTRACHET: 1                         Max.   :8516   Max.   :66897  
##  (Other)             :25                                                       
##    x_centroid     y_centroid       z_moyen        superficie     population   
##  Min.   :8288   Min.   :66477   Min.   :222.0   Min.   :  90   Min.   : 0.20  
##  1st Qu.:8350   1st Qu.:66576   1st Qu.:254.0   1st Qu.: 687   1st Qu.: 0.35  
##  Median :8450   Median :66661   Median :293.0   Median : 891   Median : 0.60  
##  Mean   :8418   Mean   :66677   Mean   :303.6   Mean   :1147   Mean   : 2.21  
##  3rd Qu.:8482   3rd Qu.:66774   3rd Qu.:349.0   3rd Qu.:1267   3rd Qu.: 1.30  
##  Max.   :8514   Max.   :66899   Max.   :464.0   Max.   :3619   Max.   :21.90  
##                                                                               
##    code_cant code_arr code_dept      nom_dept  code_reg     nom_region
##  05     :8   1:23     21:31     COTE-D'OR:31   26:31    BOURGOGNE:31  
##  24     :7   2: 8                                                     
##  15     :6                                                            
##  23     :5                                                            
##  06     :2                                                            
##  38     :1                                                            
##  (Other):2                                                            
##             Cote             geometry 
##  Côte D'or    : 3   POLYGON      :31  
##  Côte de Beane:16   epsg:2154    : 0  
##  Côte de Nuit :12   +proj=lcc ...: 0  
##                                       
##                                       
##                                       
##

Visualisation des données:

Avec ggplot, utilisation de geom_sf pour les objets spatiaux.

ggplot(data = vignobles) +
  geom_sf(aes(fill = log(PrixMoy)), lwd=0.01, colour="white")  +
  scale_fill_viridis_c(option = "viridis")+
  labs(title = "Niveaux de qualité physique des vignobles")

Pour des étiquettes réactives, on utilise geom_text_repel:

ggplot(data = communes) +
  geom_sf(aes(fill = Cote), lwd=0.01, colour="white") +
  geom_text_repel(data = communes[communes$population > 1,],
                        aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Appartenance communale aux côtes bourguignonnes")

Création de variables niveau 2 à partir d’aggrégation de données niveau 1

Info sur la distribution des prix, des surfaces, des qualités physiques des vignobles, etc.

Vignobles_par_communes <- as_tibble(vignobles) %>%
  group_by(LIBCOM) %>%
  summarise(aire_vignoble = sum(SURFACE),
            MoyPond_b = sum(SCORE_b * SURFACE) / sum(SURFACE),
            N_Parcelles = n(),
            Prix_moyen = mean(PrixMoy,na.rm=T),
            Prix_median = median(PrixMoy,na.rm=T),
            )

Info sur les appellations d’origine contrôlée:

Vignobles_par_communes_par_AOC <- as_tibble(vignobles) %>%
  group_by(LIBCOM, NIVEAU) %>%
  summarise(N_Parcelles = n()
  ) %>%  
  pivot_wider(names_from = NIVEAU, values_from = N_Parcelles) %>%
  replace_na(list(`Bourgogne` = 0,
                  `Grand cru`=0,
                  `Premier cru`=0,
                  `Village`=0,
                  `Coteaux b.`=0))

## `summarise()` has grouped output by 'LIBCOM'. You can override using the
## `.groups` argument.

Jointure et proportions d’AOC

communes_final <- left_join(Vignobles_par_communes,
Vignobles_par_communes_par_AOC, by="LIBCOM") %>%
  left_join(communes,., by=c("nom_comm"="LIBCOM")) %>%
  mutate(share_Bourgogne = `Bourgogne` / N_Parcelles * 100,
         share_GdCru = `Grand cru` / N_Parcelles * 100,
         share_PmCru = `Premier cru` / N_Parcelles * 100,
         share_Village = `Village` / N_Parcelles * 100,
         share_Coteau = `Coteaux b.` / N_Parcelles * 100,
         share_PmGrCru = (`Premier cru`  + `Grand cru`) / N_Parcelles * 100,
         share_PmGrCruVill = (`Premier cru`  + `Grand cru` + `Village`) / N_Parcelles * 100)


head(communes_final)

# Prix moyen des vins estimés par commune
ggplot(data = communes_final) +
  geom_sf(aes(fill = Prix_moyen), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="red") +
  geom_text_repel(data = communes_final[communes_final$Prix_moyen > 100 ,],
                  aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Prix moyen par commune")

# Qualité moyenne des parcelles par commune
ggplot(data = communes_final) +
  geom_sf(aes(fill = MoyPond_b), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="goldenrod4") + 
  geom_text_repel(data = communes_final[communes_final$MoyPond_b > 80
                                      | communes_final$MoyPond_b < 65 ,], 
                  aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Qualités moyennes des vignobles par commune")

# Proportion de parcelles en grand cru par commune
ggplot(data = communes_final) +
  geom_sf(aes(fill = share_GdCru), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="goldenrod") +
  geom_text_repel(data = communes_final[communes_final$share_GdCru > 8 ,],
                  aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Part de grands crus par commune")

# Proportion de parcelles en premier ou grand cru par commune
ggplot(data = communes_final) +
  geom_sf(aes(fill = share_PmGrCru), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="gold") +
  geom_text_repel(data = communes_final[communes_final$share_PmGrCru > 30 ,],
                  aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Part de premiers et grands crus par commune")

Rapatriation des données communes au niveau vignoble

communes_final_tib <- as_tibble(communes_final) %>%
  select(-geometry)

vignobles_final <- left_join(vignobles,communes_final_tib, 
                              by=c("LIBCOM"= "nom_comm"))

summary(vignobles_final)

##                              Concat                     LIBCOM    
##  1AUXEY-DURESSESCORLOIN         :   1   MEURSAULT          : 158  
##  1AUXEY-DURESSESDERRIERE LE BOIS:   1   GEVREY-CHAMBERTIN  : 128  
##  1AUXEY-DURESSESEN CHAIGNUT     :   1   BEAUNE             : 126  
##  1AUXEY-DURESSESEN FRIVOLE      :   1   NUITS-SAINT-GEORGES: 121  
##  1AUXEY-DURESSESEN MORVIN       :   1   SAINT-AUBIN        : 121  
##  1AUXEY-DURESSESEN POCHENOT     :   1   (Other)            :1727  
##  (Other)                        :2385   NA's               :  10  
##             NOM               NIVEAU        SURFACE          SCORE_b     
##  LE VILLAGE   :  26   Bourgogne  : 624   Min.   :0.0000   Min.   :18.40  
##  LES CRAIS    :  10   Coteaux b. : 319   1st Qu.:0.1300   1st Qu.:67.90  
##  LES CRAS     :   8   Grand cru  :  37   Median :0.3100   Median :73.20  
##  SOUS LA VELLE:   8   Premier cru: 410   Mean   :0.4773   Mean   :71.89  
##  LES PERRIERES:   7   Village    :1001   3rd Qu.:0.6100   3rd Qu.:76.60  
##  (Other)      :2331                      Max.   :7.2000   Max.   :94.22  
##  NA's         :   1                                                      
##             Source        PrixMoy             gid          id_geofla    
##  guide_hachette: 416   Min.   :    5.0   Min.   :  563   Min.   :  563  
##  interpolation :1942   1st Qu.:   23.0   1st Qu.: 6720   1st Qu.: 6720  
##  wine_search   :  33   Median :   36.0   Median :14392   Median :14392  
##                        Mean   :  154.4   Mean   :10484   Mean   :10484  
##                        3rd Qu.:   87.0   3rd Qu.:14597   3rd Qu.:14597  
##                        Max.   :25960.0   Max.   :14749   Max.   :14749  
##                        NA's   :4         NA's   :10      NA's   :10     
##    code_comm      insee_com                 statut       x_chf_lieu  
##  412    : 158   21412  : 158   Chef-lieu canton: 257   Min.   :8291  
##  295    : 128   21295  : 128   Commune simple  :1998   1st Qu.:8342  
##  054    : 126   21054  : 126   Sous-préfecture : 126   Median :8396  
##  464    : 121   21464  : 121   NA's            :  10   Mean   :8405  
##  541    : 121   21541  : 121                           3rd Qu.:8480  
##  (Other):1727   (Other):1727                           Max.   :8516  
##  NA's   :  10   NA's   :  10                           NA's   :10    
##    y_chf_lieu      x_centroid     y_centroid       z_moyen        superficie  
##  Min.   :66472   Min.   :8288   Min.   :66477   Min.   :222.0   Min.   :  90  
##  1st Qu.:66555   1st Qu.:8338   1st Qu.:66553   1st Qu.:256.0   1st Qu.: 765  
##  Median :66641   Median :8395   Median :66642   Median :300.0   Median :1006  
##  Mean   :66654   Mean   :8404   Mean   :66655   Mean   :306.2   Mean   :1353  
##  3rd Qu.:66758   3rd Qu.:8477   3rd Qu.:66751   3rd Qu.:355.0   3rd Qu.:1962  
##  Max.   :66897   Max.   :8514   Max.   :66899   Max.   :464.0   Max.   :3619  
##  NA's   :10      NA's   :10     NA's   :10      NA's   :10      NA's   :10    
##    population       code_cant   code_arr    code_dept        nom_dept   
##  Min.   : 0.200   05     :688   1   :1795   21  :2381   COTE-D'OR:2381  
##  1st Qu.: 0.400   15     :515   2   : 586   NA's:  10   NA's     :  10  
##  Median : 0.700   23     :457   NA's:  10                               
##  Mean   : 2.385   24     :374                                           
##  3rd Qu.: 1.600   06     :150                                           
##  Max.   :21.900   (Other):197                                           
##  NA's   :10       NA's   : 10                                           
##  code_reg        nom_region              Cote      aire_vignoble   
##  26  :2381   BOURGOGNE:2381   Côte D'or    : 152   Min.   :  5.25  
##  NA's:  10   NA's     :  10   Côte de Beane:1421   1st Qu.: 23.39  
##                               Côte de Nuit : 808   Median : 44.57  
##                               NA's         :  10   Mean   : 45.51  
##                                                    3rd Qu.: 58.18  
##                                                    Max.   :101.94  
##                                                    NA's   :10      
##    MoyPond_b      N_Parcelles       Prix_moyen       Prix_median    
##  Min.   :58.19   Min.   :  8.00   Min.   :  16.55   Min.   :  17.0  
##  1st Qu.:69.98   1st Qu.: 63.00   1st Qu.:  26.11   1st Qu.:  26.0  
##  Median :71.39   Median : 88.00   Median :  35.48   Median :  34.5  
##  Mean   :71.07   Mean   : 93.89   Mean   : 154.16   Mean   : 112.9  
##  3rd Qu.:73.59   3rd Qu.:121.00   3rd Qu.: 140.79   3rd Qu.:  84.0  
##  Max.   :87.55   Max.   :158.00   Max.   :2151.34   Max.   :1447.0  
##  NA's   :10      NA's   :10       NA's   :10        NA's   :10      
##    Bourgogne       Grand cru      Premier cru       Village     
##  Min.   : 2.00   Min.   :0.000   Min.   : 0.00   Min.   : 2.00  
##  1st Qu.:12.00   1st Qu.:0.000   1st Qu.: 8.00   1st Qu.:27.00  
##  Median :21.00   Median :0.000   Median :17.00   Median :39.00  
##  Mean   :24.49   Mean   :1.351   Mean   :16.81   Mean   :38.15  
##  3rd Qu.:34.00   3rd Qu.:2.000   3rd Qu.:26.00   3rd Qu.:48.00  
##  Max.   :55.00   Max.   :9.000   Max.   :42.00   Max.   :62.00  
##  NA's   :10      NA's   :10      NA's   :10      NA's   :10     
##    Coteaux b.    share_Bourgogne  share_GdCru      share_PmCru    
##  Min.   : 0.00   Min.   :11.11   Min.   : 0.000   Min.   : 0.000  
##  1st Qu.: 5.00   1st Qu.:19.05   1st Qu.: 0.000   1st Qu.: 7.207  
##  Median :10.00   Median :26.98   Median : 0.000   Median :19.672  
##  Mean   :13.09   Mean   :26.08   Mean   : 1.554   Mean   :17.136  
##  3rd Qu.:16.00   3rd Qu.:31.40   3rd Qu.: 2.500   3rd Qu.:24.324  
##  Max.   :41.00   Max.   :75.00   Max.   :11.111   Max.   :44.444  
##  NA's   :10      NA's   :10      NA's   :10       NA's   :10      
##  share_Village    share_Coteau    share_PmGrCru    share_PmGrCruVill
##  Min.   :11.54   Min.   : 0.000   Min.   : 0.000   Min.   :25.00    
##  1st Qu.:32.77   1st Qu.: 4.688   1st Qu.: 7.207   1st Qu.:46.84    
##  Median :39.24   Median :11.570   Median :23.529   Median :63.22    
##  Mean   :41.83   Mean   :13.398   Mean   :18.690   Mean   :60.52    
##  3rd Qu.:48.98   3rd Qu.:17.857   3rd Qu.:27.473   3rd Qu.:73.75    
##  Max.   :76.54   Max.   :44.444   Max.   :55.556   Max.   :86.49    
##  NA's   :10      NA's   :10       NA's   :10       NA's   :10       
##           geometry   
##  MULTIPOLYGON :2391  
##  epsg:2154    :   0  
##  +proj=lcc ...:   0  
##                      
##                      
##                      
##

Modélisation multiniveaux

variable a expliquer: le prix moyen du 🍷

summary(vignobles_final$PrixMoy)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     5.0    23.0    36.0   154.4    87.0 25960.0       4

ggplot(data=vignobles_final, aes(x=PrixMoy)) +
  geom_histogram(fill="coral1")+
  labs(title = "Distribution des prix du vin")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 4 rows containing non-finite values (`stat_bin()`).

# Passage en log pour une distribution plus normale
vignobles_final$LogPrix <- log(vignobles_final$PrixMoy)

ggplot(data=vignobles_final, aes(x=LogPrix)) +
  geom_histogram(fill="coral3") +
  geom_vline(aes(xintercept=mean(LogPrix)),   
             color="black", linetype="dashed", size=1)+
  labs(title = "Distribution des log de prix moyens")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 4 rows containing non-finite values (`stat_bin()`).

## Warning: Removed 2391 rows containing missing values (`geom_vline()`).

variables explicatives niveau 1:

#Score qualité geographique (pente, geologie, pedologie, precipitations):
ggplot(data=vignobles_final, aes(x=SCORE_b)) +
  geom_histogram(fill="brown1")+
  labs(title = "Distribution des niveaux de qualité")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Surface:
ggplot(data=vignobles_final, aes(x=SURFACE)) +
  geom_histogram(fill="brown2")+
  labs(title = "Distribution des surfaces de vignoble")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Niveau d'aoc:
positions <- c("Coteaux b.", "Bourgogne", "Village",
               "Premier cru", "Grand cru")
ggplot(data=vignobles_final, aes(x=NIVEAU)) +
  geom_bar(fill="brown3") +
  scale_x_discrete(limits = positions)+
  labs(title = "Distribution des niveaux d'AOC")

#Source info prix:
ggplot(data=vignobles_final, aes(x=Source)) +
   geom_bar(fill="brown4") +
  labs(title = "Distribution des sources de données par vignoble")

variables explicatives niveau 2:

#Population de la commune
ggplot(data=communes_final, aes(x=population)) +
  geom_histogram(fill="aquamarine1", bins = 5)+
  labs(title = "Distribution des populations par commune")

#Part de grand et premier cru de la commune
ggplot(data=communes_final, aes(x=share_PmGrCru)) +
  geom_histogram(fill="aquamarine3", bins = 5)+
  labs(title = "Distribution des AOC premier et grand crus par commune")

#Cote de nuit ou cote de beaune:
pos_cote <- c("Côte de Nuit", "Côte de Beane", "Côte D'or")
ggplot(data=communes_final, aes(x=Cote)) +
  geom_bar(fill="aquamarine4") +
  scale_x_discrete(limits = pos_cote)+
  labs(title = "Distribution des communes par côte")

#Qualite moyenne de la commune
ggplot(data=communes_final, aes(x=MoyPond_b)) +
  geom_histogram(fill="darkslategrey", bins = 5) +
  labs(title = "Distribution des qualités moyennes")

Centrage-réduction-sélection des variables

vignobles_final_cr <- vignobles_final %>%
  mutate(LogPrix = scale(as.numeric(LogPrix)), 
         Surface = scale(as.numeric(SURFACE)), 
         Quality = scale(as.numeric(SCORE_b)), 
         AOC = fct_reorder(NIVEAU, SCORE_b),
         Cote = fct_reorder(Cote, SCORE_b),
         MeanQuality = scale(as.numeric(MoyPond_b)), 
         ShareGdCru = scale(as.numeric(share_GdCru))) %>%
  select(LogPrix, Surface, Quality, AOC,
         MeanQuality, ShareGdCru,
         LIBCOM, Source, Cote, Concat) %>%
  na.omit


summary_stat <- data.frame()

Modèle à un seul niveau

L1 <- lm(LogPrix~Quality+Surface+AOC+Source, data=vignobles_final_cr)

summary(L1)

## 
## Call:
## lm(formula = LogPrix ~ Quality + Surface + AOC + Source, data = vignobles_final_cr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0696 -0.6275 -0.2024  0.4016  3.9311 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.71893    0.09121  -7.882 4.88e-15 ***
## Quality             -0.01648    0.03248  -0.507 0.611999    
## Surface             -0.07580    0.01951  -3.885 0.000105 ***
## AOCBourgogne        -0.15699    0.07522  -2.087 0.036995 *  
## AOCVillage          -0.20275    0.08539  -2.374 0.017655 *  
## AOCPremier cru      -0.05133    0.11279  -0.455 0.649102    
## AOCGrand cru         2.33775    0.46750   5.001 6.14e-07 ***
## Sourceinterpolation  0.98998    0.05650  17.520  < 2e-16 ***
## Sourcewine_search    0.96382    0.46841   2.058 0.039735 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8659 on 2368 degrees of freedom
## Multiple R-squared:  0.2521, Adjusted R-squared:  0.2496 
## F-statistic:  99.8 on 8 and 2368 DF,  p-value: < 2.2e-16

Comme prévu, les vignobles de grandes surfaces ont des prix en moyenne moins élevés (“ce qui est rare est cher”), et les AOC “Grands Crus” sont significativement plus chers que les AOC Côteaux Bourguignons.
En revanche, la qualité physique des vignobles est non-significative dans le modèle des prix à un niveau, et les AOC “Bourgogne” et “Village” sont significativement moins chers que les AOC Côteaux Bourguignons.

Modèle nul à 2 niveaux

mnnull <- lmer(LogPrix ~ 1 + (1 | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mnnull)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: LogPrix ~ 1 + (1 | LIBCOM)
##    Data: vignobles_final_cr
## 
##      AIC      BIC   logLik deviance df.resid 
##   4487.1   4504.4  -2240.5   4481.1     2374 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4605 -0.3560  0.0610  0.4361  4.9734 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  LIBCOM   (Intercept) 0.8638   0.9294  
##  Residual             0.3611   0.6009  
## Number of obs: 2377, groups:  LIBCOM, 31
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)   0.1520     0.1677 30.9052   0.907    0.372

print(ranova(mnnull))

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LogPrix ~ (1 | LIBCOM)
##              npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>          3 -2240.5 4487.1                         
## (1 | LIBCOM)    2 -3371.2 6746.5 2261.4  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

vrn2 <- summary(mnnull)$varcor[1][[1]][[1]]
vrn1 <- as.data.frame(summary(mnnull)$varcor)[2,]$vcov 
icc <- vrn2 / (vrn1 + vrn2)
icc

## [1] 0.7051989

La variance résiduelle du niveau commune est de vrn2 = 0.864.
La variance résiduelle du niveau vignoble est de vrn1 = 0.361.
Donc la correlation intraclasse ICC est de vrn2 / (vrn1 + vrn2) = 0.705.

Ou plus simplement:

performance::icc(mnnull)[[1]]

## [1] 0.7051989

Dans ce modèle vide à 2 niveaux, 29.5% de la variance est prise en charge par l’appartenance au niveau 2.

# Pour résumer:
summary_stat["MNNull","N"] <- dim(mnnull@frame)[1] 
summary_stat["MNNull","LogLik" ] <- summary(mnnull)$AIC[3][[1]] 
summary_stat["MNNull", "AIC" ] <-summary(mnnull)$AIC[1][[1]]
summary_stat["MNNull", "BIC" ] <- summary(mnnull)$AIC[2][[1]]
summary_stat["MNNull", "Deviance" ] <- summary(mnnull)$AIC[4][[1]]
summary_stat["MNNull", "Total_Var" ] <- vrn2 + vrn1
summary_stat["MNNull", "ICC" ] <- performance::icc(mnnull)[[1]]
summary_stat["MNNull", "InterClass" ] <- vrn1 / (vrn2 + vrn1)


summary_stat

Analyser les résidus de niveau 2

Extraction des résidus avec ranef, puis jointure et cartographie.

ranef <- ranef(mnnull)$LIBCOM %>%
  mutate(LIBCOM = as.character(rownames(.)),
         residuCom = as.numeric(`(Intercept)`))

mnnull_ranef <- left_join(communes_final, ranef, by=c("nom_comm"="LIBCOM"))

ggplot(data = mnnull_ranef) +
  geom_sf(aes(fill = residuCom), lwd=0.01, colour="white")  +
  scale_fill_gradient2(low = "#00008B", mid = "white", high = "#800020") +
  labs(title = "Résidus de niveau 2") +
  geom_text_repel(data = mnnull_ranef[abs(mnnull_ranef$residuCom) > 1,],
                        aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")

À comparer avec la carte originale des prix moyens par commune = différence de centrage-réduction.

ggplot(data = communes_final) +
  geom_sf(aes(fill = Prix_moyen), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="red") +
  labs(title = "Prix moyen du vin par commune") +
  geom_text_repel(data = communes_final[communes_final$Prix_moyen > 500 ,],
                        aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")

Modèle multiniveau avec variables de niveau 1 uniquement

mn1 <- lmer(LogPrix ~ Quality + Surface + Source + AOC +
              (1 | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mn1)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: LogPrix ~ Quality + Surface + Source + AOC + (1 | LIBCOM)
##    Data: vignobles_final_cr
## 
##      AIC      BIC   logLik deviance df.resid 
##   3038.4   3101.9  -1508.2   3016.4     2366 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4448 -0.5237 -0.0687  0.4612  4.6352 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  LIBCOM   (Intercept) 0.7839   0.8854  
##  Residual             0.1937   0.4401  
## Number of obs: 2377, groups:  LIBCOM, 31
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -0.50786    0.16599   36.31713  -3.060 0.004149 ** 
## Quality               -0.04720    0.01704 2347.41662  -2.770 0.005644 ** 
## Surface                0.04862    0.01028 2347.53697   4.728 2.40e-06 ***
## Sourceinterpolation    0.89024    0.02919 2346.53970  30.498  < 2e-16 ***
## Sourcewine_search      1.90564    0.24399 2347.88348   7.810 8.53e-15 ***
## AOCBourgogne          -0.07872    0.03881 2346.73517  -2.028 0.042646 *  
## AOCVillage            -0.17123    0.04409 2346.68467  -3.884 0.000106 ***
## AOCPremier cru        -0.07690    0.05793 2346.38460  -1.327 0.184518    
## AOCGrand cru           0.46195    0.24339 2347.85577   1.898 0.057822 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Qualty Surfac Srcntr Srcwn_ AOCBrg AOCVll AOCPrc
## Quality      0.152                                                 
## Surface     -0.010  0.229                                          
## Sorcntrpltn -0.168  0.072  0.212                                   
## Sorcwn_srch -0.006  0.060  0.191  0.094                            
## AOCBourgogn -0.197 -0.584 -0.073 -0.028 -0.028                     
## AOCVillage  -0.229 -0.720 -0.086  0.087 -0.022  0.817              
## AOCPremircr -0.237 -0.745 -0.096  0.234 -0.014  0.743  0.848       
## AOCGrandcru -0.067 -0.292 -0.258  0.002 -0.908  0.232  0.268  0.271

La relation entre prix et surface s’est inversée: les vignobles de grandes surfaces ont des prix en moyenne plus élevés.
Contre tout intuition, la qualité physique des vignobles est négativement (et significativement) associée aux prix dans ce modèle.
L’ajout des variables de niveau 1 a permis de réduire la variance totale de 20.2%.
La part de variance liée à l’appartenance communale est maintenant de 19.8%.

summary_stat["MNVar1Niv","N"] <- dim(mn1@frame)[1] 
summary_stat["MNVar1Niv","LogLik" ] <- summary(mn1)$AIC[3][[1]] 
summary_stat["MNVar1Niv", "AIC" ] <-summary(mn1)$AIC[1][[1]]
summary_stat["MNVar1Niv", "BIC" ] <- summary(mn1)$AIC[2][[1]]
summary_stat["MNVar1Niv", "Deviance" ] <- summary(mn1)$AIC[4][[1]]
summary_stat["MNVar1Niv", "Total_Var" ] <-
  summary(mn1)$varcor[1][[1]][[1]] + as.data.frame(summary(mn1)$varcor)[2,]$vcov 
summary_stat["MNVar1Niv", "ICC" ] <- performance::icc(mn1)[[1]]
summary_stat["MNVar1Niv", "InterClass" ] <- 1 -  performance::icc(mn1)[[1]]
summary_stat["MNVar1Niv", "PseudoR2" ] <- 
  (1 - (summary_stat["MNVar1Niv", "Total_Var" ]  / summary_stat["MNNull", "Total_Var" ] )) 
summary_stat

Résidus de niveau 1

mn1_residuals <- cbind(vignobles_final_cr, residuals(mn1))
colnames(mn1_residuals)[length(colnames(vignobles_final_cr))] <- "residusN1"

summary(mn1_residuals$residusN1)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.39613 -0.23049 -0.03025  0.00000  0.20294  2.03984

ggplot(data = mn1_residuals) +
  geom_sf(aes(fill = residusN1), lwd=0.01, colour="white")  +
  scale_fill_gradient2(low = "darkgreen", mid = "white", high = "darkorange") +
  labs(title = "Résidus de niveau 1") + 
  theme_minimal()

Résidus de niveau 2

ranef <- ranef(mn1)$LIBCOM %>%
  mutate(LIBCOM = as.character(rownames(.)),
         residuCom = as.numeric(`(Intercept)`))

mn1_ranef <- left_join(communes_final, ranef, by=c("nom_comm"="LIBCOM"))

ggplot(data = mn1_ranef) +
  geom_sf(aes(fill = residuCom), lwd=0.01, colour="white")  +
  scale_fill_gradient2(low = "#00008B", mid = "white", high = "#800020") +
  labs(title = "Résidus de niveau 2") +
  geom_text_repel(data = mn1_ranef[abs(mn1_ranef$residuCom) > 1,],
                        aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")

Modèle multiniveau avec variables de niveau 1 et 2

mn2 <- lmer(LogPrix ~ Quality + Surface + Source + AOC +
               MeanQuality + ShareGdCru + Cote +
              (1 | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mn2)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## LogPrix ~ Quality + Surface + Source + AOC + MeanQuality + ShareGdCru +  
##     Cote + (1 | LIBCOM)
##    Data: vignobles_final_cr
## 
##      AIC      BIC   logLik deviance df.resid 
##   2999.6   3086.2  -1484.8   2969.6     2362 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4420 -0.5218 -0.0649  0.4570  4.6376 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  LIBCOM   (Intercept) 0.1684   0.4104  
##  Residual             0.1937   0.4401  
## Number of obs: 2377, groups:  LIBCOM, 31
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -1.05896    0.11393   42.28250  -9.295 8.91e-12 ***
## Quality               -0.04788    0.01704 2346.30069  -2.809 0.005004 ** 
## Surface                0.04860    0.01028 2350.17563   4.727 2.41e-06 ***
## Sourceinterpolation    0.89129    0.02919 2347.54433  30.537  < 2e-16 ***
## Sourcewine_search      1.88415    0.24384 2353.35431   7.727 1.62e-14 ***
## AOCBourgogne          -0.07774    0.03880 2347.94052  -2.003 0.045241 *  
## AOCVillage            -0.17095    0.04408 2347.70878  -3.878 0.000108 ***
## AOCPremier cru        -0.07587    0.05793 2347.03637  -1.310 0.190419    
## AOCGrand cru           0.47943    0.24325 2352.82532   1.971 0.048845 *  
## MeanQuality           -0.04275    0.08039   30.57464  -0.532 0.598741    
## ShareGdCru             0.36251    0.08666   30.27123   4.183 0.000227 ***
## CoteCôte D'or          0.62786    0.26975   31.72140   2.328 0.026477 *  
## CoteCôte de Nuit       1.14890    0.17176   30.09307   6.689 2.04e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

La relation entre prix et surface est toujours inversée: les vignobles de grandes surfaces ont des prix en moyenne significativement plus élevés.
La qualité physique des vignobles est toujours négativement (et significativement) associée aux prix dans ce modèle.
L’AOC “Grand Cru” est positivement (et significativement) associée aux prix du vin dans ce modèle, mais moins fortement que la part de grands crus présents dans la commune (ShareGdCru).
Par rapport à la Côte de Beaune, les vins de la Côte de Nuit se vendent très significativement plus chers.
L’ajout des variables de niveaux 1 et 2 a permis de réduire la variance totale de 70.4%.
La variance interclasse est maintenant de 53.5%.

summary_stat["MNVar2Niv","N"] <- dim(mn2@frame)[1] 
summary_stat["MNVar2Niv","LogLik" ] <- summary(mn2)$AIC[3][[1]] 
summary_stat["MNVar2Niv", "AIC" ] <-summary(mn2)$AIC[1][[1]]
summary_stat["MNVar2Niv", "BIC" ] <- summary(mn2)$AIC[2][[1]]
summary_stat["MNVar2Niv", "Deviance" ] <- summary(mn2)$AIC[4][[1]]
summary_stat["MNVar2Niv", "Total_Var" ] <-
  summary(mn2)$varcor[1][[1]][[1]] + as.data.frame(summary(mn2)$varcor)[2,]$vcov 
summary_stat["MNVar2Niv", "ICC" ] <- performance::icc(mn2)[[1]]
summary_stat["MNVar2Niv", "InterClass" ] <- 1 -  performance::icc(mn2)[[1]]
summary_stat["MNVar2Niv", "PseudoR2" ] <- 
  (1 - (summary_stat["MNVar2Niv", "Total_Var" ]  / summary_stat["MNNull", "Total_Var" ] )) 
summary_stat